Monte Carlo simulations of Rb2MnF4, a classical Heisenberg antiferromagnet in 

two-dimensions with dipolar interaction 
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We study the phase diagram of a quasi-two dimensional magnetic system Rb2MnF4 with Monte 
Carlo simulations of a classical Heisenberg spin Hamiltonian which includes the dipolar interactions 
between Mn'^+ spins. Our simulations reveal an Ising-like antiferromagnetic phase at low magnetic 
fields and an XY phase at high magnetic fields. The boundary between Ising and XY phases is 
analyzed with a recently proposed finite size scaling technique and found to be consistent with a 
bicritical point at T = 0. We discuss the computational techniques used to handle the weak dipolar 
interaction and the difference between our phase diagram and the experimental results. 

PACS numbers: 68.35.Rh 75.30.Kz 75.10.Hk 75.40.Mg 
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I. INTRODUCTION 



The phase diagram of anisotropic Heisenberg antifer- 
romagnets has been studied with renormalization group 
(RG) methodsii^i^iii^ and Monte Carlo simulations,^i^i^ 
In three dimensions, RG calculations for 4 — e dimen- 
sions and Monte Carlo simulations have found an Ising- 
like antiferromagnetic (AF) phase at low magnetic fields 
and an XY phase at high fields, separated by a 1st order 
spin-flop transition line. The spin-flop transition line ter- 
minates at a bicritical point (BCP), where it meets the 
phase boundary between the XY phase and the param- 
agnetic (PM) phase, and the AF-PM phase boundary. In 
two-dimensions, due to the Mermin- Wagner theorem,— a 
BCP with 0(3) symmetry has to be at zero-temperature, 
which was confirmed by RG calculations in 2 -f e dimen- 
sions for the anisotropic non-linear cr-modeli^i^ The XY- 
PM phase boundary and AF-PM phase boundary are ex- 
ponentially close to each other while the PM phase sand- 
wiched in between narrows as exp(— 47r/T). On the other 
hand, the continuum field theory of this model contains 
an infinite number of relevant perturbations beyond the 
anisotropic nonlinear cr-model. Thus, it is also valid to 
argue that the multicritical point may not be 0(3) sym- 
metric and occurs at a finite temperature.^^ One would 
look for numerical evidence that distinguishes different 
scenarios. However, Monte Carlo simulations have been 
unable to trace the phase boundaries of the XY and AF 
phases to sufficiently low temperatures, due to the ex- 
ponentially large correlation length. Recently, a novel 
finite size scaling analysis was used to interpret the data 
from Monte Carlo simulationsi^ It was found that the 
apparent spin-flop transition line was actually consistent 
with a zero temperature BCP. An additional continuous 
degeneracy in the ground state at the spin-flop field has 
also been recently discoveredji^ The ground state actu- 
ally bears some similarities to a tetracritical phase; thus 
it was argued that the "hidden bicritical point" might be 
relabeled as the "hidden tetracritical point." 
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FIG. 1: (color online) The unit cell of Rb2MnF4 and the 
schematic phase diagram. If the BCP is at T = 0, the dashed 
line actually represents two very close second order phase 
boundaries. If the BCP is at a finite temperature, the dashed 
line represents a single first order phase transition. The the- 
oretical XY phase is found to have transverse AF order in 
neutron scattering experiments. 



In real materials, an ideal two-dimensional Heisen- 
berg spin system has not been found, since in a three- 
dimensional system, the interactions between spins can 
never be completely restricted to two dimensions. Never- 
theless, Rb2MnF4 is a very good quasi-two-dimensional 
Heisenberg antiferromagnet. In this layered compound, 
Mn^'*' ions with spin-5/2 reside on (001) planes, as shown 
in Fig. m Adjacent planes are widely separated by Rb^ 
ions, so that the exchange interactions between magnetic 
ions in different planes are negligible. The antiferro- 
magnetic order parameter has been accurately measured 
with neutron scattering experiments fl2. and analyzed 
with spin- wave theory The theoretical model with only 
nearest neighbor exchanges and a staggered magnetic 
field accounts for the experimental data very well. In 
the right hand portion of Fig. [T] we show a schematic 
phase diagram that summarizes the prevailing theoret- 
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ical alternatives and experimental data for Rb2MnF4. 
On the other hand, the large magnetic moment of Mn^^ 
ions makes it possible to model the spins with classi- 
cal vectors. Therefore, it is an excellent system to test 
theoretical predictions for two-dimensional Heisenberg 
spin systems, given that the effective anisotropy due to 
the dipolar interaction is accounted for Obviously, the 
dipolar interaction plays an important role in this sys- 
tem, as it provides the effective anisotropy that stabilizes 
the low-field AF phase and could mediate a dimensional 
crossover from two dimensions to three dimensions in the 
real material. With the in-plane isotropic exchange inter- 
action and the dipolar interaction, the Neel temperature 
at zero-field was calculated by Monte Carlo simulations 
to be 39.7±0.1 slightly higher than the experimental 
value 38.5±1.0 Kiiiiii Following the previous research, 
we performed extensive Monte Carlo simulations in both 
zero and non-zero magnetic fields to construct the full 
phase diagram and compare it with the experiments.^* 
We hope to see our model reproduce the "apparent" BCP 
at approximately T = 30K, as seen in the experiments. 
To determine the phase diagram in the thermodynamic 
limit, we used different finite size scaling analyses for dif- 
ferent phase boundaries. In particular, the "apparent" 
spin flop transition has to be examined with the novel 
finite size scaling method developed in Ref. [5, and it is 
actually found to be consistent with a zero temperature 
BCP. 

The Hamiltonian of our model reads 

n = - jsis+i) ^ s,.s,-^ ^r^^r/^f 

- ^SgnBh.- Sj, (1) 

i 

where S — 5/2, are three dimensional unit vec- 
tors, J = 0.6544meV, the dipolar interaction constant - 
U = 0.214727meVA3, the Lande .g-factor 5 = 2, the ex- 
ternal magnetic field h is fixed in the z-direction, and the 
summation over is over all nearest neighbor pairs. 
The dipolar interaction tensor V is given by: 

T^tf-i^r'^^-^-'M^^'- (2) 

The Mn^"*" ions are located on a body centered tetrago- 
nal lattice, with in-plane lattice constant a = 4.2A, and 
c-axis lattice constant c — 13.77A. However, it is known 
that the dipolar interaction between two tetragonal sub- 
lattices nearly vanishes due to the geometric arrangement 
of the moments J^!^ Therefore, besides a few simulations 
with two sublattices performed to check the validity of 
this assumption, we included only one sublattice in most 
of our simulations, which allowed us to simplify the dipo- 
lar summation and to run simulations for larger systems. 
Because the inter-layer interaction is weak, we have in- 
cluded up to four layers of spins in our simulations, with 
open boundary condition in the z direction. Each layer 
is a square lattice with lattice constant equal to a and 
the distance between adjacent layers equal to c. 



The Hamiltonian Eq. ([T]) is an approximation of the ac- 
tual quantum mechanical Hamiltonian, where spin opera- 
to rs have be en replaced with classical vector spins SSi or 
y^S{S +1)S,. Here some ambiguities arise as to whether 
S or y/S{S + 1) should be used. For the dipolar term, 
we assume that the magnetic field generated by a spin 
is a dipole field of a magnetic moment gSfiB, a-nd the 
dipolar interaction energy of a second spin with moment 
gSfiB in this field is clearly proportional to S^. This 
approximation guarantees that the total dipolar energy 
of a ferromagnetic configuration agrees with macroscopic 
classical magnetostatics of bulk materials. The exchange 
term is more ambiguous. One can argue that S'(S' -1-1) fol- 
lows from the quantum mechanical origin of the exchange 
interaction. After all, the appropriate constant should 
reproduce the correct spin wave spectrum or the criti- 
cal temperature within acceptable error bars. There is 
no guarantee that both of them can be accurately repro- 
duced with the same classical approximation. In general, 
by adopting the classical approximation to spins, one ad- 
mits an error possibly of order 1/S in some quantities. To 
justify our choice in Eq. ([1]), we first found that the criti- 
cal temperature at zero field of Eq. ^ was quite close to 
the experimental value, then we turned on the magnetic 
fields to explore the full phase diagram. It is unlikely that 
the entire experimental phase diagram would be repro- 
duced exactly including the spin-flop field. However, our 
Monte Carlo simulations should exhibit the same critical 
behavior as the real material, given that they are in the 
same universality class. In particular, we want to test if 
there is a "real" BCP at a finite temperature due to the 
long-range nature of the dipolar interaction. 

This paper is organized as the following: In Sec. [Hi 
we briefly review the simulation techniques used in this 
research, especially those designed to handle long-range, 
but very weak, dipolar interaction; in Sec. lIIIi we present 
the results from simulations performed near each phase 
boundary; in Sec. IIVI we discuss the results and give our 
conclusions. 



II. MONTE CARLO METHODS 

A. Dipole summation 

Direct evaluation of the dipolar energy in Eq. ([T]) 
should be avoided because the computational cost of di- 
rect evaluation scales as O(iV^) where N is the number 
of spins, and the periodic boundary condition needs to be 
satisfied. In our simulations we have as many as 8 x 10* 
spins and need to evaluate the dipolar energy repeatedly. 
Therefore, a fast algorithm for dipolar interaction is re- 
quired. We used the Benson and Mills algorithm^^ which 
employs the fast Fourier transformation of the spins to re- 
duce the computational cost to 0{N In N). After Fourier 
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transform, the dipolar sum in Eq. ([T]) can be written as 



E ^:f.(q)5;:(q)^f,(-q), 



(3) 



where n and n' label the different layers of the system, 
q is the in-plane wave vector, and 1?"^ / (q) is the Fourier 
transform of 2?^''. This expression is less costly to eval- 
uate than the Eq. ([2]), since the double summation over 
all the spins is replaced by a single summation over the 
wave vectors, and i^^f/(q) are constants which can be 
calculated quickly in the initialization stage of the simu- 
lation. Explicit expressions for I?"^ / (q) were first derived 
in Ref. i21i . and were reproduced in Ref. [l^ with more de- 
tail and clarity. 



B. Monte Carlo updating scheme and histogram 
reweighting 

In Monte Carlo simulations of magnetic spin systems, 
cluster algorithms offer the benefit of reduced correlation 
times. In Ref. .16i, the Wolff cluster algorithm^"^ was used 
to generate new spin configurations based on the isotropic 
exchange term in the Hamiltonian. Although the Wolff 
algorithm is rejection-free by itself, the new configura- 
tion then has to be accepted or rejected with a Metropo- 
lis algorithm according to its dipolar and Zeeman energy. 
The changes in the dipolar energy and Zeeman energy are 
roughly proportional to the size of the cluster generated 
by the Wolff algorithm. When these changes are larger 
than ksT^ the number of rejections rapidly increases, 
leading to substantially lower efhciency. This problem 
occurs when the magnetic field is typically several Tesla 
in our simulations. On the other hand, in the paramag- 
netic phase or one of the ordered phases, the cluster size 
is small, the change in dipolar energy is also small. It, 
thus, becomes redundant to evaluate the dipolar energy 
after every small change in the spin configuration. 

Since there are no rejection free algorithms for the 
dipolar interaction, and the dipolar energy only con- 
tributes a fraction of about 0.1 per cent to the total en- 
ergy in our simulations, one of our strategies to handle 
the dipolar interaction is to accumulate a series of sin- 
gle spin flips before evaluating the dipolar energy, then 
accept or reject this series of flips as a whole with the 
Metropolis algorithm depending on the change of the 
dipolar energy. The number of single spin flips for each 
Metropolis step can be adjusted in the simulation so 
that the average acceptance ratio is about 0.5, at which 
the Metropolis algorithm is most efhcient. We used the 
rejection-free heat-bath algorith m^'^'^^i^^ to perform sin- 
gle spin flips, which handles both the isotropic exchange 
and Zeeman terms in the Hamiltonian on the same foot- 
ing. 

Although fast Fourier transform significantly reduces 
the computational cost of dipolar interaction, this part 
is still the bottle-neck of the simulation. Therefore, we 



want to further reduce the number of dipolar energy eval- 
uations. To this end, we separate a short-range dipolar 
interaction from the full dipolar interaction. The short- 
range part can be defined with an cutoff in distance. In 
our simulations, we have included the up to fifth near- 
est in-plane neighbor of each spin, and the spins directly 
above or below it in the adjacent layer of the same sub- 
lattice, to form the short range dipolar interaction. This 
short-range dipolar interaction can be handled with the 
heat-bath algorithm on the same footing with the ex- 
change and the Zeeman term. The extra cost of eval- 
uating local fields produced by the additional 22 neigh- 
boring spins is insignificant. With this modification in 
single spin updates, the Metropohs algorithm should be 
performed with respect to the change in the long-range 
dipolar interaction, i.e., the difference between the total 
dipolar energy and the short-range dipolar energy. Since 
this long range dipolar energy is typically a small frac- 
tion (about 1 per cent) of the total dipolar energy, it 
is justified to accumulate many single spin flips before 
refreshing the total dipolar energy. 

We have found that the long-range dipolar energy in 
our simulations is usually a fraction of about 0.001 per 
cent of the total energy, which is actually comparable 
to ksT. This allows us to further simplify the above 
algorithm by removing the Metropolis step in the sim- 
ulation, while we simply calculate and record the full 
dipolar energy for each configuration whose energies and 
magnetizations are stored for histogram reweighting. In 
the end, we get a Markov chain of configurations from 
the simulation generated with a modified Hamiltonian 
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hort ; 



(4) 



where the the first two terms are the exchange and Zee- 
man terms in Eq. ([T]), and the last term is the short- 
range dipolar interaction. For those configurations se- 
lected for computing thermodynamic averages, we cal- 
culate and record Ti! , Tishortj their full dipolar energy 
^^dipoic, staggered magnetization of each layer 



M 



t _ 



(5) 



where L is the size of each layer and the index I is the 
layer index, and the average magnetization per spin in 
the z direction 
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L'N, 



(6) 



where TV; is the number of layers in the system. As we 
have observed that the interlayer coupling due to the 
dipolar interaction is very weak, we define the total stag- 
gered magnetization as 



N, 



1/2 



(7) 
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Similarly, the Ising-like AF order parameter is defined as 

1/2 



Ml 



(8) 



and the XY order parameter is defined as 



1/2 



(9) 



Note that we have ignored the factor SgfiB in the defini- 
tions of various magnetizations so that they are normal- 
ized to 1 in the antiferromagnetic configuration. Addi- 
tionally, the fourth order Binder cumulant for a quantity 
Q is defined as 
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(10) 



where (...) represents the ensemble average. 

The thermodynamic averages with respect to Ti.' at a 
temperature and a magnetic field slightly different from 
the simulation can be obtained with the conventional his- 
togram reweighting technique.— To calculate the ther- 
modynamic average with respect to the original Hamil- 
tonian, the weight for each sample should be modified 
to 

exp I ~ ^ ^y, ~ Sgfj,BMz{h' - h) + 7Yiong]| 



X exp 



n' 

kRTi ' 



(11) 



where Hiong — 'Tidipoic — '^shoit, T and h are the temper- 
ature and field at which the simulation was performed, 
while T' and hi are the temperature and field at which 
the histogram reweighting is done. 

The performance of this perturbative reweighting 
scheme is valid only when Tiiong is smaller or compa- 
rable to the thermal energy fc^T. For large system sizes, 
it has the same problem as the conventional histogram 
reweighting methods, i.e., the overlap of two ensembles 
defined by H, and li! decreases exponentially, leading to a 
very low efficiency. In fact, since both Tidipoio and Hshort 
are extensive quantities, we expect their difference Tiiong 
to scale as NsL? ■ Therefore, it will exceed any given ksT 
with a sufficiently large system size. For those large sys- 
tems, the above simulation scheme have to be modified 
to increase the overlap between the two ensembles de- 
fined by Ti' and Ti. Fortunately, even for our largest size 
L = 196, the long-range dipolar energy for a double layer 
system at about T = 20K and /i = 6T is mostly positive 
around 4meV, and is mostly distributed between /cbT 
and AksT. Therefore, the perturbative reweighting tech- 
nique serves to increase the weight on those configura- 
tions with lower dipolar energy, which are usually associ- 
ated with larger Ising order parameter. One might argue 
that the long-range dipolar interaction could be ignored 



since it is extremely small. Actually our simulations show 
that for the AF-PM and XY-PM phase boundaries, the 
long-range dipolar interaction is indeed negligible, but for 
the "apparent" AF-XY phase boundary its effect can be 
observed. With the perturbative reweighting technique, 
we gain knowledge of both Hamiltonians, with or without 
long-range dipolar interaction, simultaneously; hence we 
can tell where in the phase diagram the long-range dipo- 
lar interaction changes the phase boundaries. 

Most of the results presented in the next section were 
calculated with the perturbative reweighting technique, 
except part of the results for the apparent spin-flop tran- 
sition in Sec. IIII C[ where a difference larger than the 
error bar is observed. For equilibration, we ran two sim- 
ulations from different initial configurations until their 
staggered magnetizations converge within statistical fluc- 
tuations. Then each simulation ran for 5 x 10® to 2 x 10^ 
Monte Carlo steps per spin to accumulate a large amount 
of data for histogram reweighting. Early results for zero 
field were compared with simulations with Metropolis re- 
jection/acceptance steps based on the full dipolar inter- 
action; no difference larger than the error bar had been 
observed. 



III. RESULTS 

A. Low-field antiferromagnetic transition 

The zero-field AF-PM phase transition was studied 
with Monte Carlo simulations in Ref. 16, where Tc (the 
Neel temperature) was determined by extrapolating the 
crossing points of the Binder cumulant. Since we have 
adopted a slightly different model and also made a num- 
ber of changes to the Monte Carlo algorithm, we re- 
peated this calculation for testing and calibration pur- 
poses. The simulations were performed for double layer 
systems with L = 64, 96, 128, 144, 196. We also calcu- 
lated the Binder cumulant and performed finite size scal- 
ing analysia^^ with Ising critical exponents to fix the Neel 
temperature. Figure [2] shows the Ising order parameter 
(total staggered magnetization in the z-direction) for dif- 
ferent sizes at temperatures close to the Neel tempera- 
ture. Although the Ising order parameter shows a strong 
size dependence in the PM phase, the Neel temperature 
can not be determined directly from it. The Binder cu- 
mulant [^(Mj) is plotted in Fig. [3] Unlike the results 
in Ref. [l^ where the crossing points of C/4 are above all 
40K, we see in Fig. [3] that all the crossing points are be- 
tween 39. 5K and 40K. The crossing points of these curves 
move up towards the universal value of the Ising univer- 
sality class (C/4 « 0.618) as the system size increases. 
This trend is more clearly revealed by curve fitting with 
smooth splines, shown in the inset of Fig. O Because 
data points for {{MIY) and ((M])^) have smaller error 
bars, we actually did a curve fitting for those two quan- 
tities first and plotted the Binder cumulant curve with 
the fitted functions. Tc can be fixed to be between 39. 5K 
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FIG. 2: (color online) Ising order parameter(staggered mag- 
netization)for double layer systems of different sizes across 
the zero-field AF-PM phase transition. Data with full dipo- 
lar interaction do not differ from those with only short range 
dipolar interaction. 
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FIG. 3: (color online) The Binder cumulant for the Ising order 
parameter across the AF-PM phase transition at zero field. 
The inset shows a smooth spline fitting of the original data. 
Crossing points in these curves approach the Ising universal 
value(« 0.618). 



and 39. 6K, where the curves for three larger sizes cross. 
These observations suggest that the critical behavior of 
this dipolar two-dimensional Heisenberg antiferromagnet 
belongs to the Ising universality class. Therefore, we per- 
formed a finite size scaling analysis to test this prediction, 
as well as to fix the Neel temperature more accurately. 
Figure m shows the finite size scaling analysis of the Ising 
order parameter, where we plot {T/Tc — 1)L^^'^ versus 
((M|)^) L'^I^/^ , with Ising critical exponents v = 1 and 
/? = 1/8 . Clearly, all the data from different sizes fall 
nicely on a single curve. The best result is achieved by 
choosing Tc = 39.56K. Obvious deviations from a single 
curve are seen if Tc changes by O.IK, therefore we be- 
lieve the error bar for Tc is less than O.IK. Although we 
have obtained a which is only slightly smaller than 
that obtained in Ref. [l^ our data for the Ising order pa- 
rameter and its Binder cumulant are noticeably different 
from those in Ref. [l^. At the same temperature, data 
presented here are smaller than those in Ref. [l^. This 
difference is actually expected because of the difference 
in the strength of the dipolar interaction. The dipolar 
term is proportional to here in Eq. ([1]), but propor- 
tional to S{S + 1) in the previous work. 

We have also performed simulations aX h — and 
5T to study the AF-PM phase transition in a finite mag- 
netic field. The antiferromagnetic phase transition has 
been observed in both cases, but the order parameter 
changes more gradually with temperature when the mag- 
netic field is turned on. Finite size scaling with Ising ex- 
ponents have been performed. Figure [5] shows the scaling 
plot of ((Mz)^) aX h ~ 3T, which has a Hghtly lower T^. 
Long-range dipolar interaction only produces negligible 
changes in these data points. The valid regime for finite 
size scaling seems to be narrower than at h = Q, because 
some deviations are clearly seen in the low-temperature 
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FIG. 4: (color online) Finite size scaling analysis of the AF- 
PM phase transition at zero field. Data points are taken from 
Fig. [2] All of them fall onto a single curve with Ising critical 
exponents. 



data points. This could be due to the shape of the phase 
boundary, which is perpendicular to the temperature axis 
at /i = by symmetry, but not so at a finite magnetic 
field. Because of this, we change both the temperature 
and the effective anisotropy when the simulation scans 
temperature at a constant magnetic field. 
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FIG. 5: (color online) Same scaling plot as Fig. |4l but for 
simulations performed at /i = 3T. The critical temperature, 
at which the best collapsing of data points is achieved, is 
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FIG. 6: (color online) Average XY order parameter across 
the XY-PM phase boundary for double layer systems with 
different sizes. 



B. Kosterlitz-Thouless transition 

When the magnetic field is above 6T, the AF-PM phase 
transition disappears. Instead, the XY order parameter 
Eq. ([9]) becomes large at low temperatures. For a two- 
dimensional anisotropic Heisenberg antiferromagnet, one 
expects to see an XY phasCf^ii,'? in which the correlation 
function decreases algebraically. Since the dipolar inter- 
action breaks the spin rotational symmetry around the z 
axis on a square lattice, one would expect the XY phase 
to be destroyed by its presence. In case of a ferromagnetic 
model, it has been shown that above a critical strength, 
the ferromagnetic dipolar XY model exhibits a ferro- 
magnetic phase instead of an XY phase Experimen- 
tally, a "transverse" phase with long-range order has been 
found. ^® However, since the XY phase is also very sen- 
sitive to small perturbations such as crystal anisotropy 
and disorder, it is not clear whether the dipolar interac- 
tion in Rb2MnF4 alone would prevent it from entering 
the XY phase. To answer this question, we performed 
simulations in constant magnetic fields h = 6.4, 6.5 and 
7T at temperatures from 27K to 38K. Figure [S] shows 
the XY order parameter measured from these simula- 
tions for double layer systems with L — 72, 96, 128, 144, 
and 196. In all these simulations, the XY order param- 
eter increases gradually with lowering temperature in a 
broad range of temperature, and it is hard to determine 
the transition temperature from Fig. [6l They also look 
very different from the results in Ref. 8, where a transi- 
tion in the XY order parameter from zero to a finite value 
is clearly visible. There are two reasons for this. First, 
the effective anisotropy induced by dipolar interaction in 
Rb2MnF4 is very weak. The dipolar energy contributes 
only about 0.1 per cent to the total energy, while in the 
anisotropic Heisenberg model studied in Ref. the 



anisotropy is about 10 per cent to 20 per cent of the total 
energy (proportional to the anisotropy constant A). Sec- 
ondly, the magnetic field at which the simulations were 
performed (6.4T to 7T) is still close to the apparent spin- 
flop transition at about 6.2T, where the system is effec- 
tively an isotropic Heisenberg model. Experimentally, 
the existence of such an effective Heisenberg model has 
been tested Near the apparent spin-flop transition, the 
system has a large correlation length, which prevents the 
true XY critical behavior from being revealed in simula- 
tions of limited sizes. This also explains why in Fig. [6] 
((Mjy)^) increases more rapidly at 7T with decreasing 
temperature than it does at 6.5T. 

Nevertheless, we can see in Fig.[6]that the XY order pa- 
rameter decreases with system size faster at higher tem- 
peratures than at lower temperatures. In the PM phase, 
one expects the size dependence to be exponential, i.e., 
{{M^yf'^ oc exp(— 2L/^); while in the XY phase, the size 
dependence is power-law, i.e., {[Mly)^) oc L^"^^ , where 
7y is a temperature dependent exponent. On the XY- 
PM phase boundary, the critical value of this exponent 
is rjc — 1/8. Therefore, we plot {{MlyY) versus L in 
Fig. [7] with log-log scale, and try to identify the crit- 
ical temperature for the Kosterlitz-Thouless transition. 
Below the dashed line in Fig. [3 the order parameter ob- 
viously decreases faster than any power-law, which would 
be straight lines in the log-log scale. Above it, the data 
points are very close to power-law, and their slopes de- 
crease with temperature. These features are consistent 
with an XY-PM phase transition. The critical temper- 
ature Tkt is roughly 34K, estimated from Fig. [71 The 
same analysis has been done for simulations at 6.5T and 
the estimated Tkt is also near 34K. 

It has been found that if the square anisotropy is 
strong, the XY model conflrms the RG prediction that 
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FIG. 7: (color online) Log-log plot of the size dependence 
of the XY order parameter. The dashed line is a power-law 
with the critical exponent 2?7c = 1/4, used to identify the 
critical temperature. For each size, the temperatures of the 
data points are 27K, 28K,. . . , 38K from top to bottom. 



a second-order phase transition with nonuniversal criti- 
cal exponents occurs<^i^ If the anisotropy is weak, two 
possibihties for the phase diagram have been found by 
Monte Carlo simulations:^^ (1) a transition from the PM 
phase directly to the ferromagnetic phase, (2) a narrow 
XY phase is sandwiched between the ferromagnetic phase 
and the PM phase. Both of these cases might appear in 
our model if we replace the ferromagnetic phase with 
an antiferromagnetic phase. However, in all simulations 
performed above h = 6.4T, at the lowest temperature 
T = 27K, we still see that the XY order parameter de- 
creases with increasing system size. No evidence for this 
phase is evident, at least for the range of lattice size that 
could be considered. Based on this observation we believe 
if a low temperature in-plane antiferromagnetic phase ex- 
ists, it does not appear in the range of temperature and 
magnetic field where our simulations have investigated. 
Another check to exclude the transition from the PM 
phase to an Ising-like antiferromagnetic phase is to do 
the finite size scaling analysis with Ising exponents for 
the XY order parameter. We have found that it is im- 
possible to collapse all the data points in Fig. [5] onto a 
single curve, no matter what critical temperature we use. 

We have also performed simulations with a single layer 
of spins, and the results agreed with those for double layer 
systems within error bars. The results without pertur- 
bative reweighting, i.e., short-range dipolar interaction 
only, also do not differ noticeably from those with full 
dipolar interaction presented in Fig. [5] and [T] Therefore, 
we conclude that our results are consistent with an XY- 
PM transition. The main effect of the dipolar interaction 
is to provide an easy axis anisotropy, but the in-plane 
square anisotropy of the dipolar interaction is not strong 
enough to destroy the XY phase in the parameter ranges 
that we have examined. 



C. The transition from AF phase to XY phase 

Having found an Ising-like AF phase at low magnetic 
fields and an XY phase at high magnetic fields, we now 
turn to the boundary between these two phases. Pre- 
cisely speaking, we want to tell if this boundary exists 
in the thermodynamic limit, and if it exists, find where 
it is connected to the XY-PM and AF-PM phase bound- 
aries. So far, we know our system is best described by 
a two-dimensional anisotropic Heisenberg antiferromag- 
net with a very weak long-range interaction of square 
symmetry. Both the anisotropy and the long-range in- 
teraction come from the dipolar interaction. If the long- 
range component of the dipolar interaction can be com- 
pletely ignored, the XY-PM phase boundary and the AF- 
PM phase boundary meet at a zero-temperature BCP, 
as predicted by RG theory^i^ and confirmed by Monte 
Carlo simulations recently^ In this case, there is no 
real phase boundary between the XY phase and the AF 
phase. However, if the long-range component of the dipo- 
lar interaction is relevant, then the other two possibilities 
might be favored, i.e., a BCP at a finite temperature or 
a tetracritical point. In experiment, the neutron scatter- 
ing data favored a finite temperature BCP, so that the 
transition from the AF phase to the "transverse" phase 
is a first order phase transitions^ Whatever brings the 
transverse phase, which is observed to have long-range 
order, can also bring the bicritical point to a finite tem- 
perature. Because both the transverse phase and the AF 
phase have discrete symmetries, the BCP is not required 
to have a continuous (rotational) symmetry. The exis- 
tence of such a bicritical point at finite temperature does 
not violate the Mermin- Wagner theorem. 

We have performed simulations at constant tempera- 
tures T = 5, 10, 20, and 30 K and calculated both the 
Ising order parameter and the XY order parameter for 
magnetic fields between 6T and 6.4T. We found that a 
transition apparently occurs at about 6.2T at all temper- 
atures, and this transition happens over a larger range 
of magnetic field at higher temperatures than it does at 
lower temperatures. It must be pointed out that the 
location of this transition is about 0.9 to 1.1 T higher 
than the spin-flop transition in the experimental phase 
diagram. The transition field also does not show a no- 
ticeable temperature dependence, while the experimental 
spin-flop line has a positive slope. However, our result is 
in agreement with previous simulations in Rcf. 16, there- 
fore we believe this difference is a result of the classical 
approximation we have adopted and also possibly some 
other weak effects, e.g., crystal field anisotropy, that we 
have not included in our simulations. 

Figure [5] shows the Ising order parameter calculated 
at r = 20K across the transition for different system 
sizes. The left panel shows the result calculated with only 
short-range dipolar interaction, and the right panel shows 
the same data reweighted with full dipolar interaction. 
The XY order parameter which becomes large in higher 
magnetic fields is shown in Fig. [S] To tell if there is 
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FIG. 8: (color online) Ising order parameter of double layer 
systems across the apparent spin-flop transition at T = 20K. 
The data reweighted with full dipolar interaction in the right 
panel shift towards large magnetic field, and have larger error 
bars. 
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FIG. 10: (color online) Finite size scaling plot of the Ising 
order parameter at T = 20K with scaling ansatz for a zero- 
temperature BCP, where x — 1 — T* lnL/(27r) 
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FIG. 9: (color online) XY order parameter of double layer 
systems across the apparent spin-flop transition at T = 20K. 
The data reweighted with full dipolar interaction in the right 
panel shift towards large magnetic field, and have larger error 
bars. 

a BCP at a finite temperature, we need to classify the 
transition we have seen in Figs. [5] and [5] using a finite size 
scaling analysis. If it turns out to be a first order phase 
transition, a BCP must exist above 20K. The finite size 
scaling for the first order phase transition was established 
in Ref. [H. For a BCP at T = 0, Ref. i showed that 
logarithmic corrections to first order finite size scaling 
would be observed. We plot the Ising order parameter 
with the scaling ansatz for the zero-temperature BCP^ in 
Fig. [ini and with the first order scaling ansatz in FigfTT] 

In Fig.[Tni we have two tunable parameters: the critical 
field he and an effective temperature T* . The logarithmic 



corrections, powers of x = 1 — T* lnL/(27r), come from 
the spin renormalization constant calculated by RG for 
an effective anisotropic non-linear a model at T* , with 
effective anisotropy vanishing at h = he- By tuning he 
and T* , we have collapsed all the data points with short- 
range dipolar interaction onto a single curve very well. 
The data with full dipolar interaction also collapse onto 
a single curve, except for a few data points with relatively 
large error bars. Especially on the low-field side of the 
figure, the quality of collapsing is good. On the other 
hand, the first order scaling plot in Fig. [TT] shows clear 
systematic deviation in the low-field data points. This 
deviation is seen in both the left panel for short-range 
dipolar interaction and the right panel for full dipolar 
interaction. The only effect of the long-range part of the 
dipolar interaction is to shift the critical field he up by 
0.03T. Although this effect is small, it is clearly out of 
the error bars of the finite size scaling analysis. It is also 
expected from the comparison of left and right panels in 
Figs. [8] and m where the transition with the full dipolar 
interaction clearly shifts to higher magnetic fields. 

The same scaling analysis applies to the XY order pa- 
rameters as well. Figure [T^ compares two finite size scal- 
ing plots for the XY order parameter at T = 20K calcu- 
lated with short-range dipolar interaction. Obviously the 
scenario of a zero-temperature BCP fits the data better 
than a first order phase transition. 

At lower temperatures, the same scaling behavior of or- 
der parameters has been observed, and the critical field 
he turns out to be nearly identical. Figure [T^ shows the 
finite size scaling plots for Ising and XY order parame- 
ter calculated at T = lOK. Since the transition at lOK 
happens within a narrower range of magnetic field, we 
have included data points reweighted at fields different 
than that of the simulation. Data points for L = 196 
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FIG. 11: (color online) Finite size scaling plot of Ising order 
parameter at T = 20K with scaling ansatz for a first order 
phase transition, to compare with Fig. 1101 
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and zero-temperature BCP. The critical field he, and effective 
temperature T* are the same as those in Fig. 1101 



close to the transition which have large error bars are 
reweighted with different magnetic fields. Nevertheless, 
most of the data points collapse nicely onto a single curve. 
For data with short-range dipolar interactions, we have 
again found he = 6.22T; while for data reweighted with 
full dipolar interaction, the scaling plots look best if we 
choose he — 6.25T. 

Therefore, our finite size scaling so far is more consis- 
tent with a zero-temperature BCP than a finite temper- 
ature BCP above 20K. Reference d also predicts finite 
size scaling relations for the susceptibility and specific 
heat, it also predicts that the Binder cumulant C/4(Mj) 
is close to, but slightly below, 0.4 at the critical field. 



FIG. 13: (color online) Finite size scaling of the Ising (left) 
and XY (right) order parameter calculated at lOK, corre- 
sponding to a zero-temperature BCP. Data shown here are 
calculated with short-range dipolar interaction for double 
layer systems, data with histogram reweighting at different 
magnetic are also shown, he = 6.22 is the same as those in 
Fig. [121 while T* = 0.1 is smaller here. 
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FIG. 14: (color online) The Binder cumulant of the Ising 
order parameter, three curves for the larger sizes cross ap- 
proximately at ft = 6.203T and U4 = 0.54. 



We have observed the finite size scaling behavior of the 
susceptibility; however we have not seen behaviors of the 
Binder cumulant and the specific heat similar to those 
presented in Ref. 9. For the Binder cumulant. Fig. [14] 
shows that the curves for three larger sizes cross approx- 
imately at ft, = 6.203T and U4 = 0.54. This value is 
still very different from the universal value for the Ising 
universality class. However, this is actually consistent 
with the theory in Ref. [§, if one notices that here we 
have two nearly independent layers of spins. If there is 
only one layer, Ref. |9| has shown that at the critical field. 



10 



the system is effectively a single spin of length ( with 
no anisotropy, where C is the spin renormalization con- 
stant. Its angular distribution is uniform, which implies 
((M])") = l/(n + 1) and the crossing value of Ui{Ml) 
is approximately 0.4. In our simulations, since we have 
more than one layer, and they are weakly coupled, we ex- 
pect the total staggered magnetization of each layer M| is 
uniformly distributed on a sphere of radius C. Due to our 
definition of A/j in Eq. ([5]) , the distribution of is not 
a uniform distribution, although Mj^^ of each layer is dis- 
tributed uniformly. Suppose the interlayer coupling can 
be completely ignored, which is a crude approximation. 
After some simple calculations, we found the probability 
distribution of s = (Mj)^/^^ for a double layer system is 

r f, 0<s< i. 

Thus, if we ignore both the longitudinal fluctuation 
of staggered magnetization and the interlayer coupling, 
the Binder cumulant at the critical field should be 1 — 
(s^)p / (3 (s^)p). A numerical evaluation of this expres- 
sion gives 0.5334, which is very close to the crossing point 
in Fig. [131 Therefore, our simulation is consistent with 
weakly coupled multiple layers of an anisotropic Heisen- 
berg antiferromagnet. 

As for the specific heat, we have not seen a peak at 
the transition in all our simulations. Figure [15] shows 
the energy per spin and specific heat per spin calculated 
for double layer systems at T = 20K with short range 
dipolar interaction. The energy drops when the mag- 
netic field is larger than the critical field. However the 
specific heat shown in the inset does not show any sign 
of a peak. Although the error bar of the specific heat, as 
one can estimate from the fluctuation of the data points, 
is about 10 per cent, a peak which is expected to be 
similar to those discovered in Ref. |^, is clearly absent. 
However, this result is actually consistent with the finite 
size scaling theory for specific heat in Ref. ^ which shows 
that the peak in specific heat should be proportional to 
{dhc/dT)^. Because the critical field of our model is al- 
most independent of the temperature, i.e., dhc/dT w 0, 
we actually do not expect to see a peak in the specific 
heat here. 



D. Discussions 

To summarize our results, we construct a phase dia- 
gram in Fig. [TB] based on our simulations and compare it 
to the experimental phase diagram from Ref. IS. Both 
our XY-PM and AF-PM phase boundaries are close to 
experimental results, the most pronounced difference is 
the spin-flop line. Rigorously speaking, our spin-flop line 
is not a single line, but the extensions of XY-PM and AF- 
PM phase boundaries which are exponentially close to 
each other and meet at a zero-temperature BCP. The ex- 
perimental XY-AF "phase boundary" is empirical. Our 
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FIG. 15: (color online) The average energy per spin for a 
double layer system at T = 20K across the apparent spin flop 
transition. The inset shows the specific heat per spin, which 
does not show a peak similar to that of first or second order 
phase transitions. 



spin-flop line is higher in magnetic field than the experi- 
mental one and has a nearly vanishing slope, but this dif- 
ference in spin-flop field is most likely to be a consequence 
of the classical approximation which omitted quantum 
fluctuations of the spins. The anisotropic Heisenberg an- 
tiferromagnet studied in Ref. ^ offers an simple case to 
qualitatively analyze this effect. A brief derivation of the 
spin-flop fleld of this model is given in the appendix. If we 
assume the length of the classical spins is y/ S{S + 1), the 
zero-temperature spin-flop field of this simple model in 
the classical case is 4 Ja/ S{S -I- 1)(1 — A^). The spin-flop 
field of the quantum mechanical Hamiltonian is found to 
be AJS\^1 — within the linear spin- wave approxima- 
tion. More accurate results can be obtained by quantum 
Monte Carlo simulations, however, the linear spin-wave 
theory has already considerably reduced the spin-flop 
field. Since this simple model and the dipolar Heisenberg 
antiferromagnet studied here have the same critical be- 
havior near the apparent spin-flop transition, one would 
also expect the quantum effects in the latter model would 
reduce the spin-flop fleld by approximately the same 
amount. Acutally, given the classical result he ~ 6.25T, 
assuming the classical model consists of spins of length 
y/S{S + 1), the reduced spin-flop transition would be 
/ic/a/1 + l/S = 5.28T, which happens to be m agree- 
ment with the experimental value. 

Above the spin-flop line, we have observed the XY 
phase, as far as our simulations have covered, while the 
experiment shows a transverse phase. Therefore, our 
Hamiltonian certainly misses some weak but important 
effects in the real material, as the intricate correlation 
of the XY phase and the spin-flop transition is sensitive 
to many perturbations. Disorder is one of them, which 
can impose a cutoff in correlation length of the system 
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FIG. 16: Comparison between our phase diagram and the ex- 
perimental results. The experimental data points from Ref. 18 
are ploted in solid squares. 



SO that the system would not approach the ideal zero- 
point BCP from the narrow PM phase. As a result, an 
apparent finite temperature BCP would be observed and 
the apparent spin-flop transition below the "BCP" looks 
like a first order transition. The disorder can come from 
both the crystal defects and slight inhomogeneity in the 
magnetic field. The experimentally observed finite tem- 
perature BCP can also be a result of crossover to three 
dimensions due to very weak exchange between layers. 

The other facter that might have contributed to a 
phase diagram different from the experimental result 
is the exchange constant. The spin- wave analysis of 
Rb2MnF4, which provided us the exchange constant J, 
were done for systems in zero magnetic field, and the 
dipolar interaction had already been simplified to a tem- 
perature dependent staggered magnetic field acting on 
Mn^+ spins.— Therefore, the exchange integral provided 
by this theory is an effective quantity that depends on 
the particular form of the Hamiltonian which has been 
assumed. As far as we know, similar calculations have 
not been done in magnetic fields close to the spin-flop 
transition. It is not guaranteed that when the full dipo- 
lar interaction is used in the Hamiltonian, instead of an 
effective staggered magnetic field, the exchange integral 
deduced from a simplified Hamiltonian is still applicable 
and can be treated as a constant independent on either 
temperature or magnetic field. 

Finally, we show some results that justify two main 
assumptions, i.e., the inclusion of only a few layers of 
Mn^+ spins, and the omission of two sublattices. Fig- 
ure [IT] shows the Ising order parameter across the ap- 
parent spin-flop transition for systems with L = 96 but 
different number of layers. With short-range dipole inter- 
action, the result seems to saturate when we have three 
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FIG. 17: (color online) Ising order parameters calculated for 
systems at T = lOK with L — 96 and different number of 
layers. The thickness dependence is weak. It is more obvious 
in the left panel where we only include short-range dipolar 
interaction, than the right panel with full dipole reweighting. 



or more layers. After reweighting with full dipolar inter- 
action, the difference between data for different number 
of layers becomes even smaller. We estimate the change 
in he due to the change in number of layers should be 
of order O.OIT. Therefore, it is justified to do simula- 
tions with only a few layers of spins. The crossover to 
a three dimensional system will only occur at very low 
temperatures. Figure ITSbhows a finite-size scaling plot of 
the apparent spin-flop transition at T = lOK calculated 
with two sublattices. The dipolar interactions between 
two sublattices were truncated to third nearest neigh- 
bors, i.e., an Mn^+ spin feels the magnetic field generated 
by totally 32 neighboring spins in the Mn^+ layer above 
and below it belonging to the other sublattice. The mag- 
netic field contributed by spins outside this truncation 
radius should be extremely small based on our experi- 
ence with the long-range dipolar interaction. Compared 
with Fig. [131 which was calculated with a single sublat- 
tice, the difference in T* and he is negligible. We have 
enough reason not to expect the interaction between two 
sublattices to reduce the apparent spin-flop field he by 
more than O.IT. The actual additional energy due to the 
inter-sublattice dipolar interaction is found to be only 
comparable to the long-range dipolar energy. 



IV. CONCLUSIONS 

In conclusion, we have tried to explain the phase dia- 
gram of Rb2MnF4 using a classical spin model with dipo- 
lar interactions. A large amount of Monte Carlo simu- 
lations have been carried out to investigate the phase 
boundaries. Among different strategies to handle the 
dipolar interaction in the simulations, we have found 
our perturbative reweighting technique to be the most 
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APPENDIX: SPIN-FLOP FIELD AT T 
ANISOTROPIC HEISENBERG 
ANTIFERROMAGNET 



OF 



As an analogy to the dipolar Heisenberg antiferromag- 
net, we consider the simple anisotropic Heisenberg anti- 
ferromagnet with the Hamiltonian: 



FIG. 18: (color online) Finite-size scaling plot for simula- 
tions at T = lOK with inter-sublattice dipolar interactions. 
The data in this figure only differ very slightly from those in 
Fig. 1131 in which the inter-sublattice dipolar interactions were 
omitted. 



suitable for very weak dipolar interactions in Rb2MnF4. 
The phase diagram inferred from our data captures the 
main features of the experimental phase diagram and the 
agreement is good at low magnetic fields. On the appar- 
ent spin-flop line, the XY and AF boundaries come so 
close together that they cannot be distinguished below 
an "effective" BCP at T w 30K. However, our data anal- 
yses support a zero temperature BCP. This conclusion 
is based on a novel finite size scaling analysis for two- 
dimensional anisotropic Heisenberg antiferromagnetsi^ If 
this multicritical point is located at very low finite tem- 
perature, as suggested by Ref. [HI. We believe its tem- 
perature must be sufficiently low, which is beyond our 
numerical accuracy. The ground state degeneracy for 
the anisotropic Heisenberg antiferromagnets as found in 
Ref. [13 may also exist in our model with dipolar inter- 
actions, which we have not yet verified. If it exists, one 
might simply rename the bicritical point as a tetracritical 
point. The zero temperature BCP is located above the 
experimental spin-flop line in the phase diagram, which 
appears to be a a line of flrst order phase transitions. We 
believe this difference from the experimental phase dia- 
gram is mainly caused by the classical approximation. 
Nevertheless, we have confirmed that the dominant ef- 
fect of the dipolar interaction in Rb2MnF4 is to provide 
an effective anisotropy, while other effects, such as in- 
plane square anisotropy and interlayer interaction, are 
extremely weak. Therefore, we would hope to obtain a 
more accurate phase diagram if we performed quantum 
Monte Carlo simulations for a simpler Hamiltonian which 
includes the effective anisotropy. 



H = J ^ [A {SfS^ + SfS^) + s:s^] • - ^ E 

(A-1) 

which is defined on a square lattice. The classical version 
of this model has been well studiedj^i^i^ where the spins 
are treated as unit vectors. The spin-fiop field at zero- 
temperature is He — 4J\/1 — A^. If we replace the spins 
with vectors of length ^/S{S+^), He is then modified 
to AJy/ S{S + - A2). This is not the only way to 
make connections to the quantum Hamiltonian. One can 
also replace J with JiS'(iS' +1), while replacing H with 
HS, which can be justified by arguing that the Zeeman 
energy of the ferromagnetic configuration takes on the 
correct macroscopic value. In this case, the spin-flop field 
is modified to He = 47(5 -I- l)vl— A^. However, in 
any case, we will show that the classical spin-flop field 
is larger than the quantum mechanical spin-flop field. 
By introducing the Holstein-Priniakoff (HP) bosons on A 
and B sublattices respectively, and keeping the quadratic 
terms, the Hamiltonian Eq. (|A.1[) can be rewritten as 

^ ^ [A5(al6]+c.c.) +(4a,-5)(5-6]6,)' 

-HY,4a + Hj2blb, (A.2) 

where a and are HP boson operators on sublattice A, 
b and on sublattice B, index i labels sites on sublattice 
A which are nearest neighbors of the sites on sublattice 
B labeled with j. After a Fourier transformation, this 
quadratic Hamiltonian turns out to be 

n = -AJS{S + 1)Na - HNa - ^k, (A.3) 

k 

where Na is the number of sites on sublattice A, and 
H.^5.(„. (A,4) 
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For simplicity, we have defined h = H/SJ and 7k = 
2cosA:2; + 2cosA:j,. The spin- wave spectrum can be ob- 
tained with the Bogohubov transformation: 



Ck — cosh ^k^k + sinh 9kb^_ 



k' 



c?k — sinh^kflk + cosh^k^-k- 



(A.5) 
(A.6) 



In order to ehminate the cross terms in the Hamiltonian, 
one sets tanh26'k = A7k/4. Apart from a constant term, 
the spin- wave part of the Hamiltonian turns out to be 



where 



uj± (k) ^ JSJW - A272 ± H. 



(A.7) 



(A., 



When H is large enough such that ui- (0) becomes neg- 
ative, the AF ground state becomes unstable since the ex- 
citations on spin- wave mode Ck=o lower the ground state 
energy. This precisely indicates the the spin-flop insta- 
bility. Therefore, the critical magnetic field is given by 



Although the above spin- wave analysis is only a crude ap- 
proximation, we see that the quantum effect lowers the 
spin-flop field by a factor oiS/{S+l) or ./S/JsTT), de- 
pending on which classical approximation one uses. The 
case with A = 2/3 and S* = 1/2 has been studied with 
quantum Monte Carlo simulations,^^ Its phase diagram 
shows the spin-flop field is at approximately h/ Jxy = 1-8, 
i.e.. He = 1.2 J in our notation here. The above spin- 
wave approximation gives He = 1.49 J, and the two clas- 
sical approximations gives He = AAIJ and He = 2.58J 
respectively. Clearly, the classical approximations over- 
estimate the spin-flop field. The difference from the real 
spin-flop field is large as we expect the quantum fiuctua- 
tion to have a strong effect for S = 1/2. For larger spins, 
such as 5' = 5/2 which is studied in this paper, the classi- 
cal approximation should work better. However, we still 
expect it to overestimate the spin-flop field by an notice- 
able amount. 



He = AJSVl - A2. 



(A.9) 
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